# Sources:
# 1. http://de.wikipedia.org/wiki/Nationalratswahlordnung
# 2. http://commons.wikimedia.org/wiki/File:Regionalwahlkreise_Oesterreich.png

######################### Load Libraries #######################################

library(foreign)

library(seatsvotes)

############################# Functions ########################################

# modified version of Linzer's swr function for use with pr systems

swrpr <- function (fit, frac, seats, sims = 1000, graph = TRUE) 
{
    starttime <- Sys.time()
    dat <- reconstruct(fit)
    dat[is.na(dat)] <- 0
    np <- ncol(dat)
    votemat <- NULL
    seatmat <- NULL
    for (s in 1:sims) {
        vsim <- sim.election(fit, dat)
        vsim[is.na(vsim)] <- 0
        partyvotes <- vsim[, 1] * vsim[, -1]
        votemat <- rbind(votemat, colSums(partyvotes)/sum(partyvotes))
        seatmat <- rbind(seatmat, pr(vsim,frac,seats))
    }
    ret <- list()
    truevotes <- dat[, 1] * dat[, -1]
    ret$truevote <- colSums(truevotes)/sum(truevotes)
    ret$trueseat <- pr(dat,frac,seats)
    # names(ret$trueseat) <- names(ret$truevote)
    swing <- NULL
    swing.lm <- NULL
    seatlist <- list()
    for (j in 1:(np - 1)) {
        smean.hi <- mean(seatmat[(votemat[, j] < (ret$truevote[j] + 
            0.012)) & (votemat[, j] > (ret$truevote[j] + 0.008)), 
            j])
        smean.lo <- mean(seatmat[(votemat[, j] > (ret$truevote[j] - 
            0.012)) & (votemat[, j] < (ret$truevote[j] - 0.008)), 
            j])
        swing <- c(swing, round((smean.hi - smean.lo)/0.02, 2))
        swing.lm <- c(swing.lm, round(coefficients(lm(seatmat[, 
            j] ~ votemat[, j]))[2], 2))
        vrange <- seq(min(floor(votemat[, j] * 100)), max(floor(votemat[, 
            j] * 100)), 1)/100
        seatlist[[j]] <- cbind(vrange, matrix(NA, nrow = length(vrange), 
            ncol = 3))
        colnames(seatlist[[j]]) <- c("vote", "mean", "lower", 
            "upper")
        for (i in 1:length(vrange)) {
            sel <- seatmat[(votemat[, j] >= (vrange[i] - 0.002)) & 
                (votemat[, j] < (vrange[i] + 0.002)), j]
            if (length(sel) < 10) {
                seatlist[[j]][i, 2:4] <- NA
            }
            else {
                seatlist[[j]][i, 2] <- mean(sel)
                seatlist[[j]][i, 3:4] <- quantile(sel, probs = c(0.025, 
                  0.975))
            }
        }
    }
    if (graph) {
        par(mar = c(4.4, 4.1, 0.5, 0.5), mfrow = c(1, 1))
        pdf("swrpr_aus.pdf")
        plot(ret$truevote[!is.na(swing)], ret$trueseat[!is.na(swing)], 
            pch = 19, cex = 1.2, xlim = c(0, 0.7), ylim = c(0, 
                0.7), xlab = "Party vote share", ylab = "Party seat share", 
            cex.lab = 1.5, xaxt = "n", yaxt = "n")
        axis(1, seq(0, 1, 0.1), seq(0, 1, 0.1), cex.axis = 1.5)
        axis(2, seq(0, 1, 0.1), seq(0, 1, 0.1), cex.axis = 1.5)
        for (j in 1:(np - 1)) {
            if (!is.na(swing[j])) {
                lines(seatlist[[j]][, 1], seatlist[[j]][, 2], 
                  lwd = 2.5)
                lines(seatlist[[j]][, 1], seatlist[[j]][, 3], 
                  lwd = 2, lty = 3)
                lines(seatlist[[j]][, 1], seatlist[[j]][, 4], 
                  lwd = 2, lty = 3)
                text(ret$truevote[j] + 0.01, ret$trueseat[j], 
                  paste(colnames(dat)[j + 1], ": ", swing[j], 
                    sep = ""), cex = 1.2, pos = 4)
            }
        }
        dev.off()
    }
    vdiff <- as.data.frame(t(t(votemat) - colMeans(votemat)))
    sdiff <- as.data.frame(t(t(seatmat) - colMeans(seatmat)))
    swing.median <- round(apply(sdiff/vdiff, 2, median), 2)
    names(seatlist) <- colnames(votemat) <- colnames(seatmat) <- names(ret$truevote)
    names(swing) <- names(swing.lm) <- names(swing.median) <- names(ret$truevote)
    swingtab <- cbind(round(100 * ret$truevote, 1), round(100 * 
        ret$trueseat, 1), swing)
    rownames(swingtab) <- colnames(dat[, -1])
    colnames(swingtab) <- c("Votes", "Seats", "Swing ratio")
    print(swingtab)
    ret$votemat <- votemat
    ret$seatmat <- seatmat
    ret$swing <- swing
    ret$swing.lm <- swing.lm
    ret$swing.median <- swing.median
    ret$vs <- seatlist
    ret$time <- Sys.time() - starttime
    return(ret)
}

# add swrpr function to Linzer's seatsvotes package
environment(swrpr) <- asNamespace("seatsvotes")

# new pr function modeling Austrian electoral rules 
pr <-
function(votes,frac,seats) {
	
	votes <- votes[,1]*votes[,-1]
	# votes <- (votes[,1]*10000)*votes[,-1]
	votes <- colSums(votes)
	votes <- as.data.frame(rbind(votes,matrix(NA,frac,length(votes))))
	rownames(votes) <- NULL
	nc <- ncol(votes)
	divisor <- seq(1+1,frac+1,by=1)
	
	for (j in 1:nc) {

		votes[-1,j] <- votes[1,j]/divisor

	}

	seatrank <- as.vector(apply(votes,2,cbind))

	names(seatrank) <- as.vector(apply(as.matrix(colnames(votes)),1,rep,frac+1))

	seatrank <- sort(seatrank,decreasing=TRUE)

	wm <- NULL

	wm <- c(wm,which(rank(-seatrank)<=seats))
	
	tb <- table(names(wm))/seats

	return(tb[names(votes)])

}

# add pr function to Linzer's seatsvotes package
environment(pr) <- asNamespace("seatsvotes")

############################ Data Setup ########################################

########################## Austria #############################################

# Documentation: aus_elections.xls

# change to directory with relevant electoral data
setwd("~/_data/electoral/aus")

# generate vector of election years
aus_yrs <- c(1945,1949,1953,1956,1959,1962,1966,1970,1971,1975,1979,1983,1986,1990,1994,1995,1999,2002,2006,2008)

# generate empty list with length equal to the number of election years
AUS.elections <- vector("list",length(aus_yrs))

# use election years as names for items of list 
names(AUS.elections) <- sprintf("%04.0f",aus_yrs)

# read csv files for election years into R and assign the data matrices to respective items of the elections list
for (i in aus_yrs) {

	name <- paste(i, "election", sep = "_")
	
	yr <- paste(i,sep ="")

	AUS.elections[[paste(yr)]] <- read.csv(paste(name,"csv", sep="."), sep=",")
	
}

# remove missing data
for (i in 9:15) {

	AUS.elections[[i]] <- AUS.elections[[i]][is.na(AUS.elections[[i]]$year)==FALSE,]
	
}

# use districts as row names; eliminate unnecessary columns; rename first column 
for (i in 1:length(AUS.elections)) {

	rownames(AUS.elections[[i]]) <- AUS.elections[[i]]$district_name
		
	AUS.elections[[i]] <- AUS.elections[[i]][-which(names(AUS.elections[[i]]) %in% c("year","district_name","X"))]
	
	
	names(AUS.elections[[i]])[1] <- "Total"	
}

############################ Analysis ##########################################

########################## Austria #############################################

################################## SM Figure 2 #################################

# set random number seed for replicability
set.seed(1234)

# set working directory for graphs; a pdf of the graph will be generated as part of the swrpr function 
setwd("~/_supportingmaterials_tabfig/")

# Steps:
# 1. Generate data frame for relevant election year.
# 2. Find contestation patterns (seatsvotes packages).
# 3. Estimate mixture model based on contestation patterns identified in Step 2 (seatsvotes package).
# 4. Verify fit of mixture model estimated in Step 3 (seatsvotes package).
# 5. Simulate seats-votes elasticities (seatsvotes package using modified functions - see above).

##########
## 2006 ##
##########

# Step 1
aus06 <- AUS.elections$'2006'
aus06 <- data.frame(Total=aus06$Total,
                   OVP=aus06$vsh_ovp,
                   FPO=aus06$vsh_fpo,
                   GRUNE=aus06$vsh_grune,
                   SPO=aus06$vsh_spo,
                   BZO=aus06$vsh_bzo)

# seats won: SPÖ: 68; GRÜNE: 21; FPÖ: 21; BZÖ: 7; ÖVP: 66

# Step 2
aus06.split <- findpatterns(aus06)

# Step 3
fit06 <- list()
fit06[[1]] <- mvnmix(dat=aus06.split[[1]],components=4,nrep=10,scatter=T)

# Step 4
show.marginals(fit06,numdraws=50000)

# Step 5
sr06 <- swrpr(fit06, 99, 183, sims=10000)

